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Proximity  control  of  modem  nano-spacecraft  often  relies  on  low  and  discrete  thrust 
engines  that  are  characterized  by  low  consumption,  and  generate  on-off  force  profiles. 
New  guidance  solutions  must  take  into  account  the  nature  of  this  type  of  orbital  engines. 
This  paper  introduces  novel  analytical  guidance  solutions  for  spacecraft  relative  motion 
considering  continuous,  on-off  thrust,  and  using  relative  orbit  elements  as  a  geometrical 
representation  of  the  dynamics.  The  solutions  provide  the  relative  state  vector  at  any 
given  time,  accommodating  any  thrust  magnitude  along  the  three  directions  of  the 
relative  frame,  as  well  as  generic  activation  times  and  durations.  Relative  orbit  elements 
geometrically  interpret  key  aspects  of  the  relative  motion,  including  for  example,  the 
relative  ellipse  size,  and  the  evolution  of  its  center  in  time.  The  new  solutions  provide  the 
guidance  designer  with  a  direct  visualization  of  the  thrust  effects  on  the  relative  motion 
geometry,  offering  new  possibilities  for  analytical  guidance  in  the  presence  of  continuous 
thrust  engines,  such  as  low  thrust  engines  on  nano-spacecraft.  The  paper  presents  the 
analytical  solutions,  and  tests  their  effectiveness  using  a  sample  thrust  profile  based  on 
input-shaping,  previously  developed  by  one  of  the  authors  using  classical  Cartesian 
coordinates.  The  use  of  relative  orbit  elements  shows  substantial  benefits  and  added 
simplicity  with  respect  to  Cartesian-based  approaches,  holding  the  promise  for  straight¬ 
forward  onboard  spacecraft  implementation.  The  software  developed  for  this  research 
will  be  available  open  source1,  to  be  used  by  spacecraft  guidance  designers  as  trajectory 
design  tool. 

©  2014  IAA.  Published  by  Elsevier  Ltd.  All  rights  reserved. 


1.  Introduction 

Small  spacecraft  flying  in  tight  formations  are  nowadays  replacing  larger  single  satellites,  due  to  their  lower  cost,  the 
reconfiguration  ability,  the  flexibility  to  substitute  malfunctioning  vehicles  without  aborting  the  mission,  and  their  inherent 
redundancy  as  multi-vehicle  systems  [1].  On  the  other  hand,  solutions  such  as  the  CubeSats2,  present  a  new  set  of  design 
challenges,  mainly  related  to  the  vehicles’  limited  size,  power,  and  computation  abilities.  Incorporating  thrusters  and 
carrying  on-board  propellant  is  extremely  difficult  on  nano-spacecraft  weighting  a  few  kilograms  [2],  and  such  thrusters 
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Nomenclature 

Q 

orbital  right  ascension  of  the  ascending 

node  (RAAN) 

a 

orbital  semi-major  axis 

(Op 

orbital  argument  of  perigee 

ae 

relative  ellipse  semi-major  axis 

^Earth 

mean  radius  of  the  Earth 

Ax,y,z,i 

amplitude  of  control  force  at  the  ith  firing  in 

ROEs 

relative  orbit  elements 

the  x,  y,  or  z  direction 

s 

Laplace  complex  variable 

P 

parametric  (phase)  angle  for  the  planar 

t 

time 

motion  in  terms  of  ROEs 

t* 

unknown  time  variable  in  input-shaping 

C 

reference  thrust  value  (it  represents  the  nom- 

profile 

inal  thrust  available  on  a  spacecraft) 

x,y,  z 

Cartesian  coordinates  in  the  LVLH  frame 

AtXtytz,i 

duration  of  the  ith  firing  in  the  x,  y ,  or  z 

xd,yd 

center  of  the  2-by-l  relative  ellipse  (part  of 

direction 

the  ROEs  variables) 

A  tw 

unknown  duration  in  input-shaping  profile 

x,  y,  z 

Laplace  transforms  of  the  Cartesian  coordi¬ 

e 

orbital  eccentricity 

nates  in  the  LVLH  frame 

r 

parametric  (phase)  angle  for  the  out  of  plane 

motion  in  terms  of  ROEs 

Further  symbols  explanation 

HCW 

Hill-Clohessy-Wiltshire 

*orb 

orbital  inclination 

Subscript  0  refers  to  initial  conditions  (at  initial  epoch 

h 

Earth’s  second  zonal  harmonic 

time  t0).  Subscript  /  refers  to  final  epoch. 

LVLH 

local  vertical  local  horizontal 

Dot  on  a  variable  represents  first  time  derivative.  Two 

P 

Earth’s  gravitational  constant 

dots,  second  time  derivative,  etc. 

n 

orbital  angular  rate 

Subscript  h  indicates  coasting  solutions  to  the  relative 

N*,y,z 

number  of  firings  in  the  x,  y,  or  z  direction 

motion  dynamics. 

V 

orbital  polar  angle 

operate  at  one  -  or  just  a  few  -  nominal  value  of  force,  i.e.  they  are  on-off  only.  As  for  the  computational  capabilities,  very 
simple  programs  must  be  designed  for  the  vehicles  to  be  autonomous.  Analytical  solutions  are  needed  for  straightforward 
online  implementation,  and  to  completely  avoid  the  need  of  onboard  numerical  iterations. 

The  relative  motion  of  spacecraft  formations  is  commonly  represented  in  a  relative  frame  using  Cartesian  coordinates. 
Relative  orbit  elements  (ROEs)  represent  a  nonlinear  transformation  from  Cartesian  coordinates  to  geometric  variables, 
giving  a  visual  and  straightforward  understanding  of  the  main  aspects  of  proximity  flight  dynamics.  Other  researchers  have 
presented  various  solutions  separating  the  oscillatory  and  drifting  motions  in  the  classical  linearized  equations  of  spacecraft 
relative  motion  [3-6],  using  linear  transformations.  These  previous  efforts  are  not  directly  and  thoroughly  addressing  the 
geometrical  problem  of  relative  motion.  In  particular,  ROEs  are  akin  to  classical  orbital  elements,  in  that  they  consist  of 
physical  lengths  and  angles  allowing  easy  visualization  of  any  relative  orbit  (a  benefit  not  provided  in  Refs.  [3-6]). 

This  paper  presents  the  general  analytical  solution  for  the  time  evolution  of  the  ROEs,  when  on/off  constant  thrust  is 
used.  These  results  are  of  particular  interest  for  missions  employing  low  thrust  engines.  The  new  solutions  also  hold  the 
potential  for  on-board  implementation.  Alternately,  given  their  analytical  nature,  they  may  serve  as  an  initial  guess  for 
numerical  optimizers  to  minimize  fuel/time,  and  enable  verification  of  various  pre-designed  thrust  profiles.  In  this  paper  the 
authors  demonstrate  the  last  feature,  by  deriving  solutions  for  orbital  planar  re-phasing  (moving  to  a  new  location  along 
track)  or  rendezvous  (moving  to  the  location  of  a  chief  satellite,  i.e.  the  origin  of  the  relative  motion  frame)  using  thrust 
profiles  based  on  input-shaping. 

Input-shaping  has  been  extensively  used  in  vibration  suppression  for  flexible  manipulators  Refs.  [7-14],  but  never  for 
orbital  control,  to  the  authors’  knowledge.  Input-shaping  is  a  convolution  technique  based  on  the  knowledge  of  a  system’s 
natural  frequencies  of  oscillation.  Given  a  feed-forward  control  signal,  designed  to  perform  a  desired  maneuver,  but  not 
taking  into  account  potential  excitation  of  undesired  oscillations,  input-shaping  consists  of  the  convolution  of  the  signal 
itself  and  a  specified  train  of  impulses,  so  that  the  system’s  resulting  behavior  presents  minimal  residual  vibrations  at  the 
end  of  the  maneuver.  The  impulses  and  their  locations  in  time  are  computed  based  on  the  frequencies  that  need  to  be 
suppressed,  i.e.  the  modes  one  wants  to  limit  in  amplitude.  The  majority  of  input-shaping  applications  falls  under  the 
category  of  flexible  structures  control,  such  as  space  manipulators  control.  It  is  important  to  underline  that  input-shaping  is 
not  intended  to  reduce  the  energy  of  a  system.  Roughly  speaking,  existing  oscillations  cannot  be  damped  with  input¬ 
shaping,  while  maneuvers  from  an  equilibrium  set  to  a  new  equilibrium  set  are  possible,  as  in  the  case  of  re-phasing 
maneuvers.  In  the  specific  context  of  spacecraft  relative  motion,  oscillations  refer  to  periodic  motion  in  the  position 
coordinates. 

Exploiting  the  new  analytical  formulas,  the  special  case  of  an  input-shaping  profile  is  presented,  and  the  analytical 
solution  for  spacecraft  planar  rendezvous  with  along-track  control  only  is  derived.  In  addition,  the  paper  demonstrates  how 
the  input-shaped  control  profile  can  be  ad-hoc  modified  to  obtain  a  final  close  relative  motion  of  desired  size  relative  to  a 
reference  satellite.  Sample  numerical  simulations  show  some  of  the  maneuvers  achieved  via  the  analytical  solutions. 
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The  intended  advancements  in  the  state  of  the  art  for  relative  motion  guidance  design  are: 

•  Use  of  ROEs  in  conjunction  with  on-off  thrust  profiles,  thus  enabling  geometrical  visualization  of  the  key  aspects  of  relative 
motion,  and  addressing  modern  engines  found  especially  in  small  satellites. 

•  Analytical  solutions  for  ROEs  time  evolution,  and  examples  of  their  use  with  a  specific  open-loop  thrust  signal  and  a 
closed-loop  application. 

•  Illustration  of  potential  future  uses  for  the  new  analytical  formulas. 

The  paper  is  organized  as  follows.  Section  2  presents  the  spacecraft  relative  motion  dynamics  in  Cartesian  coordinates  and  its 
nonlinear  transformation  in  ROEs.  Section  3  is  dedicated  to  the  derivation  of  the  general  analytical  equations  for  the  ROEs 
evolution  in  time  when  on-off  thrust  is  used.  Section  4  shows  the  example  where  an  input-shaping-based  along-track  thrust 
profile  is  applied  to  the  new  analytical  equations,  to  derive  close  form  guidance  solutions  for  re-phasing  maneuvers.  Section  5 
illustrates  the  guidance  obtained  in  the  previous  section  with  numerical  simulations.  The  same  section  also  presents  one  closed- 
loop  example  where  the  guidance  is  computed  iteratively  when  used  in  a  more  realistic  nonlinear  simulation  environment. 
Section  6  draws  the  conclusions  and  suggests  future  applications  for  the  new  analytical  solutions.  The  software  developed  in 
Matlab®  and  Simulink®  for  this  investigation  will  be  made  available  open  source  (link  in  Ref.  [15]),  for  interested  researchers  and 
guidance  designers. 

2.  Satellite  relative  dynamics 

Consider  two  satellites  orbiting  in  close  proximity  to  each  other.  For  this  analysis,  one  will  be  referred  to  as  the  reference 
satellite,  or  “chief,”  and  the  other  as  the  “deputy.”  For  the  methods  presented  here,  it  is  assumed  that  the  only  force  acting 
on  each  satellite  is  that  of  a  point  mass  gravitational  field,  the  chief  is  in  a  circular  orbit,  and  the  distance  between  the 
satellites  is  small  compared  to  their  orbital  radius.  These  assumptions  yield  the  following  linear  time-invariant  differential 
equations  [16,17]: 

x-2ny-3n2x  =  0 
y+2nx  =  0 

z+n2z=  0  (1) 

These  are  known  as  the  Hill-Clohessy-Wiltshire  (HCW)  equations  and  are  written  in  the  local-vertical,  local-horizontal 
(LVLH)  coordinate  frame,  whose  origin  is  at  the  chief  satellite.  In  these  equations,  x  is  the  component  of  the  deputy’s  position 
vector  relative  to  the  chief  in  the  radial  direction  positive  away  from  the  Earth,  y  is  the  along-track  component  positive  along 
the  velocity  vector  of  the  chief,  and  z  is  the  cross-track  component  perpendicular  to  the  orbital  plane  of  the  chief,  n  is  the 
mean  motion  of  the  chief.  The  LVLH  frame  is  depicted  in  Fig.  1. 

The  solution  to  Eq.  (1)  is: 


z  =  —  sin  (nt)+z0  cos  (nt) 


z  =  —  sin  (nt)+z0  cos 


n 


x  =  x0  cos  (nt)  +  (3nx0  +  2y0)  sin  (nt) 
y  =  -  2x0  sin  (nt)  +  (6 nx0 + 4y0)  cos  (nt)  -  (6nx0  +  3y0) 
z  =  z0  cos  (nt)-nzo  sin  (nt) 


(2) 


x  (radial) 


Fig.  1.  Depiction  of  LVLH  Frame. 
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where  x0,  y0,  etc,  are  conditions  at  some  epoch  time  t0,  and  t  is  the  time  since  t0.  Consider  the  following  change  of 
coordinates  from  x,y,z,x,y,z  [18]: 

ae  =  2'l/(£)2  +  (3x+2i)2  xd=4x+2l 
yd=y-2*j3  =  atan2(x,  3  nx + 2  y) 

zmax  =  J  +z2y  =  atan2  (nz,  z)  -  atan2(x,  3nx+ 2y)  (3) 


where  ae ,  xd ,  yd,  A  zmax,  and  y  are  the  ROEs.  The  inverse  of  this  transformation  is 
x  =  ^js-  cos  p+xd  x  =  ^n  sin/? 

3 

y  =  ae  sin  /3+yd  y  =  aen  cos  (3--nxd 

z  =  zmax  sin  0 +/?)  z  =  zmaxn  cos  (y  - +/?)  (4) 

It  has  been  shown  in  Ref.  [18]  how  the  ROEs  evolve  with  time: 


flg  —  Cleo  xd  —  xd0 

yd=ydo-|rWdot=ydo-|n^t  P  =  Po  +  nt  (5) 

^max  —  Zmax  0  Y  —  Yo 

These  equations  are  analogous  to  Eq.  (2)  for  x,y,z,x,y,z  in  that  they  express  the  ROEs  values  at  any  given  time  as  a 
function  of  their  initial  (epoch)  values  and  the  time  since  epoch. 

The  parameterization  of  Eq.  (4)  reveals  that  the  relative  motion  of  the  deputy  with  respect  to  the  chief  in  the  x-y  plane  is 
a  superposition  of  periodic  motion  in  x  and  y,  with  period  equal  to  that  of  the  chiefs  orbit,  and  secular  motion  in  y. 
Essentially,  this  is  an  elliptical  path  that  is  drifting  in  the  y-direction  at  a  rate  of  -3/2 nxd.  The  instantaneous  center  of  the 
ellipse  is  ( xd,yd ).  It  has  a  semi-major  axis  of  length  ae  in  the  along-track  direction  and  semi-minor  axis  of  length  ae\2  in  the 
radial  direction.  (3  is  a  parametric  angle  (i.e.  phase  angle)  indicating  the  location  of  the  deputy  satellite  in  its  trajectory,  with 
/3=0  corresponding  to  the  perigee  location  (the  “bottom”  of  the  ellipse).  The  relative  motion  in  x  and  y,  if  the  elliptical  path 
were  “frozen”  at  a  point  in  time,  is  depicted  in  Fig.  2.  Although  the  ellipse  is  actually  drifting,  it  has  been  frozen  in  order  to 
conveniently  label  the  ROEs.  The  z-component  of  the  relative  motion,  according  to  the  HCW  model,  is  purely  sinusoidal 
and  independent  of  x  and  y.  This  motion  is  a  simple  harmonic  oscillator  with  amplitude  zmax  and  phase  angle  y+/?.  The 
deputy  intersects  the  chiefs  orbit  plane  at  y+/?= 0  and  k,  and  reaches  zmax  and  -zmax  at  y+/?=;r/2  and  3n/2,  respectively. 
Thus,  y  represents  the  phase  difference  between  the  x  and  y  motion  and  the  z  motion.  Fig.  3  depicts  a  typical  3-D  relative 
trajectory,  with  zmax  and  y  labeled.  (NOTE:  Because  (3  and  y  are  angular  representations  of  time  -  similar  to  mean  anomaly  - 
they  are  labeled  in  Figs.  2  and  3  as  /?*and  y*  which  are  the  physical  interpretations  of  these  angles.) 


orbital 

velocity 

direction 


Fig.  2.  Planar  projection  of  relative  motion  trajectory  with  relative  orbit  elements  labeled. 
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relative 


Fig.  3.  Depiction  of  out-of-plane  relative  motion  with  relative  orbit  elements  labeled. 


Fig.  4.  Generic  example  of  on-off  continuous  thrust  profile. 


3.  Analytical  solutions  with  continuous  on-off  thrust 

This  section  presents  the  steps  to  derive  the  closed-form  solutions  for  the  time  evolution  of  the  ROEs  when  a  generic 
on-off,  continuous  thrust  profile  is  assumed  in  each  direction  of  the  LVLH  reference  frame.  AXytZ  i  indicates  the  magnitude  of 
the  ith  firing  in  the  x,  y,  or  z  direction.  Atx,y?z>l  is  the  corresponding  time  duration  of  the  firing,  while  is  the  coasting  (off) 
time  duration  between  the  (i  - 1  )th  and  the  ith  firing.  tF  is  the  final  time  (see  Fig.  4).  Note  that,  if  the  first  firing  in  a  particular 
direction  begins  at  t=  0,  then  t/i  in  that  direction  is  defined  to  be  0. 

Because  the  dynamics  we  started  from  are  linear  (Eq.  (1)),  the  superposition  principle  can  be  applied  to  find  the  state  at 
the  final  time.  In  particular,  the  final  state  can  be  written  as  the  sum  of  the  value  at  the  final  time  when  coasting  from  the 
initial  condition,  plus  each  of  the  final  values  obtained  by  starting  at  zero  initial  conditions,  coasting  for  a  duration  equal  to 

-1  (t/.  +  Atj) +t/.,  applying  the  generic  ith  thrust  for  its  given  duration,  and  then  coasting  for  a  duration  equal  to 

tF  ~'Zj  =  i  (Afj +  */,  )• In  the  previous  expressions  the  subscript  indicating  the  direction  of  the  firing  was  removed,  indicating 
its  validity  for  any  axis.  For  each  of  the  x,  y,  and  z  components,  the  ROEs  offer  a  simple  solution,  since  coasting  from  a  set  of 
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initial  conditions  is  represented  by  the  equations: 

Qe  =  Qeo 
*d=*dO 

yd=ydo~inxd0tF 

P  =  Po  +  ntF  (b) 

Zmax  —  ^max  0 
Y  =  Y0 

Finding  the  final  state  after  firing  from  zero  initial  conditions,  and  then  coasting,  requires  the  combination  of  Cartesian 
coordinates  to  find  the  state  right  after  firing,  then  conversion  to  ROEs,  and  finally  coasting,  using  the  same  form  as  in 
Eq.  (6).  Eq.  (11 )  gives  the  values  of  the  Cartesian  relative  states  after  a  generic  single  firing  of  duration  At  and  coasting  period 
t  beforehand,  with  components  in  the  x,  y,  and  z  directions.  They  can  be  derived  using  Laplace  transform  on  the  system  in 
Eq.  (1),  when  applying  control  accelerations  Ax>y<z<i.  If  s  is  the  Laplace  complex  variable,  we  start  from: 

s2X(s)  -  sx0  -  x0  -  2 n  [s Y(s)  -  sy0]  -  3 n2X(s)  =  ^ 
s2Y(s) -sy0 -y0+2n[sX(s)  -x0]  =  A 

s2Z(s)-sz0-z0  +  n2Z(s)  =  y  (7) 

Solving  Eq.  (7)  in  the  Laplace  domain  we  obtain: 


X(s)- 

 1 

s2-3n2  -2ns 

^+sx0+x0-2nsy0 

_Y(5)_ 

—  s2(s2  +  n2) 

2  ns  s2 

"^+syo+yo  +  2nXo 

Z(s)(s2  +  n2)  =f+sz0 +z0 


which  further  simplifies  into: 


X(s)=Xh(s)  + 

Y(s)  =  Yh(s)  + 

A 


Ax 


-+ 


2nAv 


s(s2  +  n2)  s2(s2  +  n2) 


4AV 


2nAx 


3AV 


Z(s)  =  - 


s(s2  +  n2)  s2(s2  +  n2) 
+Z/,(s) 


"s(s2  +  n2) 

and  finally,  converting  back  in  the  time  domain: 
x  =  xh+U  1-  cos  (nt)]+ 2^ [t 


(9) 


sin  (nt) 


y=Yh+4^[1-  cos  (nt)]- 2-^  t- 
nz  n 

z  =  zh+p1-  cos  (nt)] 

x  =  xh+y  sin(nt)  +  2y[l  -  cos  (nt)] 

y=y/i+4^  sin  (nt)  -  2y[l  -  cos(nt)]-3/\yt 

z  =  zh-\ — -  sin  (nt) 
n 


-jV 


(10) 


where  xh,  yh ,  and  represent  the  solution  of  the  HCW  equations  for  unforced  motion  (i.e.  Eq.  (2)).  To  apply  the  super¬ 
position  principle  described  earlier,  we  only  need  to  retain  the  portion  of  Eq.  (10)  generated  by  control  accelerations, 
i.e.  we  consider  null  initial  conditions.  This  provides  Eq.  (11). 


x+=£§[l-  cos(nAt)]  +  2^[t-^9] 


_  /Am  . 


m  c  ( ri  A  _ "1 


sin  (nAt) 


A  f2 
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x+  sin(nAt)  +  2^[l  -  cos(nAt)] 

y+  =4^  sin(nAt)-2^[l  -  cos(nAt)]-3AyAt 


z+=^sin(nAt)  (11) 

Eq.  (3)  is  then  used  to  convert  the  Cartesian  relative  states  (Eq.  (11))  into  ROEs,  and  the  ROEs  are  propagated  for  the 
coasting  period  according  to  Eq.  (6).  This  is  repeated  for  each  single  firing,  with  Nx,  Ny,  Nz  indicating  the  total  number  of 
firings  along  each  axis.  By  adding  together  all  the  states  obtained  as  described  above,  the  following  closed-form  solutions  for 
the  ROEs  subject  to  generic  thrust  profiles  are  obtained: 


ae(tF )  =  2 


1 


(ae0  /2)  sin  (/J0  +  ntF) 

-  t  \/-(2/n4)A^(-l+  cos (riAtXl))  sin  (  -0£  -n(tF-  £  (AtXj  +  tfXj) 


2=1 

N, 


j=  1 


+2  2  \J ~ (2/n4)^.(— 1  +  cos (nAtyi))  sin  +n^tF-  ’Z(Atyj  +  tfy.) 
(Qe0  /2)  cos  (0O + ntF) 

-  2  \j — (2/n4)AX( ( —  1  +  COS (riAtx,))  cos  ( -Aj-nUf-  £  (AtXi  +  tfXj) 


i  =  1 

N, 


J  =  1 


+2  2  \/-(2/n4y^.(-l  +  cos (n4ty,))  cos  +n(  tf-  2  (Atyj+tfy.) 


j  =  i 


/?+  =  atan2((AXi /n)  sin (nAt*,),  -(A*,/n)(l-  cos(nAtXi))) 
/?+  =  atan2((2Ay,/n)(l  -  cos  (nAty.)),  (2Ay./n)  sin  (nAty.)) 


/2\ 

Xd(tf)  =  Xdo  +  ^-J  2 

yd(tF)=ydo- (J)nxdotf- ZAXiAtXl- g) 

-3  S  AyiAtyf  ^tF  —  (^fy,  +  t/y,) ^ 


0(tF)  =  a  tan  2 


(l/2n) 


ae„n2  sin(/?0  +  ntf)-  \  \ 

2  2A<, ^/(2  —  2  cos  (n4tXj))  sin  (  -0+  -ntF+n  £  (Atx.  +  t/x.)  j  - 


i  =  1 
N, 


J=1 


2  44y.y(2-2  cos  (nAty.))  sin  (  -/?+  -ntF+n  £  (Atyi+tfy.) 


2=1 

/  Qe0n2  cos(/?0  +  ntF)  + 


J  =  1 


(2/2n) 


2  \J(2-2  cos  (nAtXi))  cos  -ntp+n  'Z  (AtXj+tfXj)  j  + 

2  4Ayiy/(2-2  cos  (nAty,))  cos  f  -0+  -ntf  +  n  2  (4tyj.  +  t/y.) 


j  =  i 


// 


/?+  =  atan2((AXi /n)  sin  (nAtX|),  —  (>lX|/n)(l  -  cos (nAtXi)) 
=  atan2((2/\y./n)(  1  -  cos  (nAtyi)), (2Ayi /n)  sin  (nAtyi)) 


(12a) 


(12b) 


(12c) 
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Zmax(£f)  — 


\ 


f  zmax0  COS  (r0 +/^0  +  +  \ 

E  (Az, /n2)y/(2-2  cos(nAtz,))  cos  (  - ¥i  -  n  (  tF  -  £  (AtZj.  +  t/z.) 


j  =  i 


^max0  sin  (ro+/?o  +  ntF)- 


2  (AZi/n2)<J (2-2  cos (nA tz,.))  sin  [  tf-  £  (AtZj+t/z.) 


j'=i 


VAj  =  atan2((^Zi/n)(l  -  cos  (nAtZi)),  (AZi/n)  sin  (nAtZi) 
/nZmaxo  sin  (r0+^0  +  nff)- 


y(tF)  =  atan2 


\ 


N; 


2=1 

NZmaxo  COS  (/q  +  A)  +  ntf)  + 
N; 


E  (Aji /n)  1^2 - 2  cos (n A tZ( )  sin  (  tF-  E  (AtZj.  +  t/Zj)  j  j, 


E  (AZj/n)^2-2  cos(nAtZ|)  cos  (  -^,-nuF-  2  (AtZj+tfz.) 


-m) 


i//i  =  atan2((^z./n)(l  -  cos  (nAtZi)),  (AZi/n)  sin (nAtZ/)) 


(12d) 


Eqs.  (12a)-(12d)  were  obtained  through  a  combination  of  symbolic  calculation  and  numerical  verification  in  Matlab®. 
Roughly  speaking,  each  firing,  i.e.  each  non  zero  phase  in  the  example  of  Fig.  4  corresponds  to  one  instance  of  Eq.  (11),  its 
transformation  into  ROEs  (Eq.  (3)),  followed  by  coasting  until  final  time  (Eq.  (6)).  The  difference  between  each  firing  is 
represented  by  its  duration  (At  in  Eq.  (11))  and  the  time  remaining  to  reach  final  time  (tF  in  Eq.  (6)  becomes 
tF  -  Ej  =  i  (Afj  +  f/j),  with  the  i  and  j  indexes  explained  earlier  in  this  section.  The  use  of  software  tools  enabled  compact 
formulation  of  the  final  analytical  solutions  (12a)-(12d),  and  the  scripts  used  to  obtain  and  validate  them  will  be  available 
open  source. 

Despite  their  complicated  appearance,  Eqs.  (12a)-(12d)  represent  a  powerful  tool  for  trajectory  design,  since  they 
are  analytical  and  because  they  heavily  simplify  for  specific  applications.  The  input-shaping  example  in  the  following 
section  shows  one  such  simplification,  and  in  general,  real  spacecraft  applications  may  reduce  the  number  of  variables  in 
(12a)-(12d),  for  example  having  only  one  value  of  thrust,  or  fixed  durations  of  the  firings,  etc. 


4.  Example  of  planar  application  of  the  roe  formulas:  Input-shaping  thrust  profile 

In  this  section  a  validation  of  some  of  Eqs.  (12a)-(12d)  is  performed.  In  particular,  one  of  the  results  previously  obtained 
by  one  of  the  authors  using  Cartesian  coordinates  Ref.  [19]  is  confirmed  by  means  of  ROEs,  obtaining  a  simpler  expression.  In 
Ref.  [19],  an  input-shaping-based,  y-only  thrust  profile  was  proved  to  be  an  effective  means  to  obtain  analytical  leader- 
follower  re-phasing  or  rendezvous  guidance,  as  well  as  equilibrium-relative-orbit  to  equilibrium-relative-orbit  guidance. 
Such  a  profile  allows  for  in-plane  control,  moving  the  center  of  the  ellipse  to  a  new  desired  location,  where  the  ellipse 
collapses  to  a  point  for  leader-follower  maneuvers.  The  thrust  profile  was  presented  in  Ref.  [19]  as  follows:  (12a) 

^  =AZi=  0  u  =  c  x  sign (ydo  -ydf),  c  >  0  tF  =  3t*  +  2 Atw 

AU, 5, 6  =  ±  4u’  tf i, 2,4,6  =  °»  AtC. ...6  —  T  (13) 

A3  4  =  ±  \ U ,  tf35  =  Atw 

Section  5  shows  the  typical  shape  of  the  input-shaping  profiles.  Representative  experiments  showing  how  input-shaping 
can  be  applied,  for  example,  to  bang-bang  control  profiles  can  be  seen  in  the  video  in  Ref.  [20].  The  profile  of  Eq.  (13) 
consists  of  known  amplitudes  for  the  firings  (c  is  a  given  control  amplitude),  while  the  A tw  and  t*  are  to  be  determined. 
ydo  and  yd  are  the  initial  and  final  (desired)  along  track  positions  of  the  relative  ellipse’s  center,  respectively.  Substituting 
Eq.  (13)  into  Eq.  (12b)  and  assuming  ydo  >yd  ,  the  following  expressions  are  obtained: 

Xd(tF)  =  xdo 

yd(tf)=yd0-(|)c(r*)2  (14) 

Eq.  (14)  leads  to  the  solution  for  t*,  given  initial  and  desired  final  values  foryd,  that  is,  initial  and  final  centers  of  the 
ellipse  of  relative  motion. 

//do-yd(tf) 

V  3  u 


(15) 
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Note  that  this  result  is  not  as  straightforward  to  find  in  Cartesian  coordinates  (Ref.  [19]),  in  which  case  there  is  also  no 
geometrical  interpretation. 

Substitution  of  the  profile  of  Eq.  (13)  in  Eq.  (12a)  does  not  lead  to  an  expression  of  comparable  simplicity.  Nevertheless, 
several  observations  can  be  made  that  provide  useful  insight  with  regards  to  the  expected  final  value  for  ae.  First  of  all,  all 
the  terms  where  thrust  along  x  appears  are  zero.  Secondly,  the  terms  not  containing  aeo  in  Eq.  (12a),  in  the  square  powers, 
represent  modifications  with  respect  to  the  initial  value  of  ae.  In  fact,  if  no  thrusting  was  present,  the  final  value  for  ae  would 
be  aeo,  as  expected.  These  observations  justify  focusing  on  only  some  of  the  resulting  terms  in  Eq.  (12a),  and  specifically  we 
here  analyze  the  following  portion,  where  the  square  power  is  omitted  for  simplicity: 

2  .Z  ^-(Jr)4(-1+  cos(nAty,)) 

sin  (p*  +  n  -  .  Z  ( Af> 'i  +  l!y, )  <16) 

After  some  algebra,  and  the  use  of  Prosthaphaeresis  formulas,  Eq.  (16)  becomes: 

/  2  sin  (J5y  +nt* +  n At w)  cos  (nt* +  n At w)+  \ 

+2  sin  (p+  +  n|t*  +  nAtw)  cos(nt*  +  nAtw)  +  ^ 

^  +4  sin  (/?+  +n|t*  +  nAtw)  cos  (n£)  ^ 

where  the  /?+  become  a  common  /?+,  given  the  nature  of  the  firings  of  same  duration  in  the  profile  of  Eq.  (13).  Eq.  (17)  still 
provides  little  information  about  what  to  expect  at  the  end  of  the  firing  sequence.  Since  t*  is  determined  in  Eq.  (15),  as  well 
as  the  /?+,  through  Eq.  (12c),  the  only  free  variable  in  Eq.  (17)  is  the  wait  time  between  the  series  of  firings  A tw.  One 
observation  to  be  made  is  that  the  term  under  the  square  root  is  never  expected  to  be  zero,  since  it  would  imply  firing  with 


i.efi- 


n4  V 


cos  |  n— 


Table  1 

Initial  Orbital  parameters  for  S/C  and  desired  trajectory  for  Leader-Follower  case,  plus  general  data  for  simulations. 


Initial  orbital  parameter 

Chief 

Deputy 

Semi-major  axis  a 

6778.1  km 

6778.1  km 

Eccentricity  e 

0 

0 

Inclination  iorb 

97.9908  deg 

97.9908  deg 

Right  ascension  of  the  ascending  node  (RAAN)  Q 

261.621  deg 

261.621  deg 

Argument  of  perigee  mp 

30  deg 

30  deg 

Polar  angle  v 

27.216  deg 

27.18  deg 

Additional  parameters  used  for  the  simulations 

^ Earth  =  6378.1363  km  ^  =  398,600.4418  km3/s2 

Fig.  5.  Example  of  ae  vs.  A tw  for  leader  follower  initial  condition.  Note :  The  above  graph  is  obtained  using  the  numerical  data  of  Table  1,  showing  the  min 
and  max  points. 
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no  duration.  For  this  reason  we  need  only  focus  on  the  parenthesis  term.  The  derivative  of  this  parenthesis  term  with  respect 
to  A tw  yields 

2  cos  +nAtw+|nt*)  +2  cos  +  nAtw+|nt*)  + 

+  2  cos  (/?+  +  2nAtw  +  2nt*^  +2  cos  +  2nAtw+^nt*^  (18) 

This  shows  that  at  the  most  four  values  for  Atw  can  represent  a  minimum/maximum  for  Eq.  (17),  within  an  orbital  period 
(0  <  Atw  <  T).  In  fact,  such  a  derivative  is  composed  of  four  cosine  functions,  all  shifted  by  different  phases. 

The  locations  of  these  minimum/maximum  points  change  from  case  to  case,  depending  on  the  values  of  /?+  and  t*. 
Despite  the  impracticability  of  solving  Eq.  (12a)  in  terms  of  Atw,  even  when  simplified  with  the  input-shaping  profile,  the 
derivative  information  allows  us  to  predict  the  type  of  function  we  should  expect,  and,  in  addition,  Eq.  (17)  clearly  shows  a 
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Fig.  6.  Rendezvous  with  chief  starting  from  an  initial  relative  point.  Top:  (1)  Atw=0.5  T,  exact  rendezvous  with  chief;  (2)  Atw=0,  obtaining  the  maximum 
ae  for  the  final  equilibrium  orbit  around  the  chief;  (3)  Atw=0.25  T,  obtaining  an  intermediate  value  of  ae  for  the  final  equilibrium  relative  orbit  around  the 
chief.  Bottom  3  plots:  control  profiles. 
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content  in  frequency,  when  using  A tw  as  independent  variable,  not  exceeding  2 n.  The  Nyquist-Shannon  sampling  theorem 
Ref.  [21]  enables  capturing  the  nature  of  the  function  representing  ae  when  input-shaping  is  applied,  and  Atwis  the 
independent  variable,  by  computing  Eq.  (12a)  only  at  At  points  spaced  by  a  l/4n  distance,  that  is,  theoretically  8^  (i.e.  25  or 
more)  points  total  in  one  orbital  period  time  frame.  A  desired  ae  value  can  be  then  interpolated  using  these  required  values 
(e.g.  using  splines),  or  more  points,  for  increased  accuracy  purposes,  posing  no  computational  issues. 

Depending  on  the  initial  conditions,  the  extrema  for  the  ae  value  can  be  four  or  less,  and  located  at  different  A tw  values 
between  0  and  the  orbital  period  T,  as  shown  later  on.  In  all  cases  there  are  special  values  of  A tw  that  zero  out  the  increase  in 
ae,  that  is,  there  are  no  oscillation  size  increases  due  to  performing  the  maneuver. 

It  should  be  noted  that  for  the  other  term  under  the  radical  in  Eq.  (12a),  an  identical  expression  can  be  found,  the  only 
difference  being  that  the  sine  function  in  Eq.  (16)  would  be  replaced  by  a  cosine  function.  Thus,  the  analysis  of  this  term 
would  be  quite  similar  to  that  above. 


5.  Sample  numerical  simulations 

In  the  following  numerical  simulations  we  assume  a  chief  satellite  located  at  the  origin  of  the  LVLH  frame,  and  that  we 
are  maneuvering  a  deputy  satellite.  The  chief  represents  the  target  trajectory  for  the  different  types  of  maneuvers  here 
presented,  i.e.  we  set  up  rendezvous  problems.  More  generally,  such  a  target  trajectory  can  be  a  virtual  satellite,  and  can  be 
located  anywhere  such  that  the  chief  and  deputy  orbital  periods  are  equal.  The  following  numerical  simulations  are 
obtained  using  the  results  presented  earlier.  One  example  of  closed-loop  control  is  also  presented  where  the  ROEs-based 
guidance  is  recomputed  when  reaching  its  final  time,  for  three  times.  This  improves  accuracy  when  the  proposed  guidance 
is  used  with  the  more  realistic  nonlinear  Keplerian  dynamics  plus  J2,  and  provides  a  proof  for  potential  flight 
implementation.  For  all  the  simulations  the  control  value  c=2x  10“ 5  m/s2  is  used,  representing  a  low-thrust  thruster. 
In  principle,  any  c  value  can  be  chosen,  representing  the  thrust  available  on  the  spacecraft. 

The  initial  orbital  parameters  of  Table  1  are  used  to  generate  the  trajectories  for  the  first  simulation,  representing  an 
initial  condition  of  leader-follower.  Note  that  the  initial  orbital  parameters  are  first  converted  to  Cartesian  position  and 
velocity  in  an  Earth  centered  inertial  frame,  then  translated  kinematically  into  the  LVLH  frame,  and  finally  forced  to  match  a 
leader  follower  initial  condition  for  the  linear  equations,  i.e.  cancelling  any  residual  relative  velocity  and  x  displacement. 


Table  2 

Initial  Orbital  parameters  for  S/C  and  desired  trajectory  for  equilibrium-to-equilibrium  case. 


Orbital  parameter 

Chief 

Deputy 

Semi-major  axis  a 

6778.1  km 

6778.1  km 

Eccentricity  e 

0 

0.0001 

Inclination  iorh 

97.9908  deg 

97.9908  deg 

Right  ascension  of  the  ascending  node  (RAAN)  Q 

261.621  deg 

261.621  deg 

Argument  of  perigee  cop 

30  deg 

30  deg 

Polar  angle  v 

27.216  deg 

27.18  deg 

Atw[s] 


Fig.  7.  Example  of  ae  vs.  A  tw  for  equilibrium  relative  orbit  initial  condition.  Note :  The  above  graph  is  obtained  using  the  numerical  data  of  Table  2,  showing 
the  min  and  max  points. 
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Fig.  8.  Rendezvous  with  chief  starting  from  an  initial  relative  equilibrium  orbit.  Top:  (1)  Atw=0.5T,  obtaining  an  intermediate  ae  (between  initial  and 
maximum  achievable)  on  final  relative  orbit;  (2)  Atw=625  s,  obtaining  the  minimum  ae  for  the  final  equilibrium  orbit  around  the  chief;  (3)  Atw=4440  s, 
obtaining  the  maximum  of  ae  for  the  final  equilibrium  relative  orbit  around  the  chief.  Center:  zoom  of  the  final  relative  orbits.  Bottom  3  plots:  control 
profiles. 
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Three  different  final  conditions  are  chosen  for  this  simulation,  one  being  exact  rendezvous  and  two  cases  where  the  final 
motion  is  a  relative  closed  orbit  around  the  chief.  For  these  cases,  the  variation  of  final  ae  as  function  of  A tw  reduces  to  a 
simple  cosine  function,  with  maximum  at  Atw=0  and  one  orbital  period,  and  no  increase  at  one-half  orbital  period  (see 
Fig.  5).  Fig.  6  shows  the  resulting  trajectories  applying  input-shaping,  as  well  as  the  control  profiles  as  dictated  by  Eqs.  (13) 
and  (15).  Note  that  in  each  case,  the  motion  is  simulated  beyond  tF  {tF  is  indicated  on  the  control  plots  for  each  case  in  Fig.  6). 
This  is  done  to  illustrate  clearly  the  final  trajectory  achieved  in  each  case. 

Table  2  introduces  a  small  eccentricity  in  the  deputy  initial  orbital  parameters,  thus  creating  an  initial  motion  which  is  a 
relative  closed  orbit  whose  center  is  offset  from  the  chief  by  the  same  amount  as  the  leader-follower  separation  in  the 
previous  cases.  Note  that  the  initial  orbital  parameters  are  first  converted  to  Cartesian  coordinates  in  an  Earth  centered 
inertial  frame,  then  translated  kinematically  into  the  LVLH  frame,  and  finally  forced  to  match  an  equilibrium  motion  initial 
condition  for  the  linear  equations,  i.e.  imposing  the  condition  y0  =  -2wx0  Ref.  [17].  For  these  scenarios,  the  final  ae  function 
is  more  complicated  than  before.  Fig.  7  indicates  that  Atw=625  s  yields  no  change  in  ae,  Atw=4440  s  yields  the  maximum 


Fig.  9.  Graphical  demonstration  of  the  number  of  points  needed  to  represent  the  ae(Atw)  function. 


y  [m] 


Fig.  10.  Example  of  closed-loop  guidance  solving  the  ROEs  input-shaping-based  solution  iteratively,  with  Atw=0.  (Top)  first  iteration,  (bottom)  and  two 
more  iterations. 
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final  value  of  ae,  and  Atw=0.5  T  yields  an  intermediate  final  value  of  ae.  These  results  are  shown  in  Fig.  8,  including  both  the 
x-y  trajectories  and  the  control  profiles.  Note  again  that  the  motion  is  simulated  beyond  tF  in  each  case  to  illustrate  the  final 
trajectory  achieved. 

All  the  maneuvers  can  be  computed  analytically,  from  Eq.  (15)  and  the  earlier  observations  on  the  function  ae(Atw). 
The  only  numerical  operation  required  to  design  such  maneuvers  consists  of  reconstructing  ae(Atw)  by  means  of  computing 
Eq.  (12a)  at  a  few  points,  and  interpolating  when  a  desired  change  in  ae  is  given,  to  solve  for  the  corresponding  A tw.  This 
provides  a  powerful  tool  to  design  guidance  trajectories  onboard  spacecraft  with  limited  computing  capabilities. 

Fig.  9  compares  three  reconstructions  of  the  ae(A£w)  function:  the  one  using  the  necessary  8 n  points,  minimally  differing 
from  the  more  accurate  line  obtained  with  a  sample  time  of  1  s.  The  third  line  shows  how  less  than  8^  points  (10  in  the 
example)  lead  to  a  poor  reconstruction  of  the  curve.  The  circles  indicate  the  (Atw,  ae )  points  required  for  the  curve 
reconstruction.  Once  those  are  stored  in  a  table,  a  desired  ae  value  leads  to  the  corresponding  Atw  by  linearly  interpolating 
between  the  two  closest  ae  points. 

Fig.  10  shows  an  example  where  the  ROEs  input-shaping-based  solution  for  equilibrium  relative  orbit  to  equilibrium 
relative  orbit  is  solved  iteratively,  to  obtain  a  closed-loop  simulation  considering  a  more  realistic  nonlinear  dynamics  for  the 
relative  motion.  The  solution  of  Eq.  (15)  with  Atw=0  is  recomputed  three  times,  at  the  end  of  each  sequence.  The  maneuver 
is  intended  to  move  the  center  of  the  already  excited  equilibrium  relative  orbit,  starting  from  the  same  initial  conditions  of 
the  simulation  presented  in  Fig.  8  and  Table  2.  The  nonlinear  relative  motion  is  simulated  using  Keplerian  dynamics  plus  J2 
for  each  satellite,  and  then  projecting  the  relative  position  and  velocity  vectors  in  the  LVLH  frame.  Fig.  10  shows  that  the  first 
iteration  achieves  a  position  magnitude  error  of  445  m  between  the  analytically  (re)generated  guidance  and  the  nonlinear 
trajectory.  This  error  actually  increases  in  the  second  iteration,  but  then  decreases  to  a  very  acceptable  56  m  in  the  third 
iterationa  decrease  in  the  error  between  the  analytically  (re)generated  guidance  and  the  nonlinear  trajectory,  starting  from 
445  m  and  ending  at  56  m  after  3  iterations.  Particularly,  the  bolded  lines  in  the  bottom  plot  of  Fig.  10  highlight  the  guidance 
trajectory  and  corresponding  nonlinear  trajectory  for  the  third  iteration,  showing  a  close  match  between  the  two  motions 
for  the  entire  duration  of  the  trajectory.  In  performing  this  simulation,  it  was  noted  that  the  accuracy  with  which  the 
nonlinear  trajectory  achieved  the  desired  final  closed  relative  orbit  was  consistent  with  the  position  magnitude  error 
described  above.  That  is,  the  error  between  the  actual  final  ROEs  achieved  by  the  nonlinear  trajectory  and  the  desired  final 
ROEs  increased  from  the  first  to  the  second  iteration,  and  decreased  from  the  second  to  the  third  iteration.  In  a  real  mission 
scenario,  mid-course  corrections  may  be  advised,  to  maintain  a  lower  tracking  error. 

More  generally,  the  new  equations  providing  the  ROEs’  time  evolution  in  analytical  form,  when  continuous,  on-off  thrust 
is  applied,  hold  the  potential  for  testing  and  designing  new  open  loop  control  sequences.  They  could  also  provide  analytical 
initial  guesses  for  numerical  optimization  of  the  guidance. 

6.  Conclusion 

This  paper  presents  the  general  analytical  solutions  for  spacecraft  relative  orbit  control,  when  on/off  continuous  thrusters 
are  used,  employing  relative  orbit  elements  instead  of  classical  Cartesian  coordinates  to  represent  the  relative  dynamics. 
Relative  orbit  elements  are  a  powerful  tool  to  visualize  geometrical  aspects  of  spacecraft  relative  motion.  A  thrust  profile 
based  on  the  input-shaping  technique  is  used  to  validate  the  obtained  formulas.  The  analytical  solutions  for  exact  re-phasing 
or  rendezvous  using  input-shaping  are  provided,  along  with  the  expressions  and  procedures  to  control  the  size  of  the  final 
relative  orbit  around  the  target  trajectory  or  chief  satellite.  Sample  numerical  simulations  show  the  type  of  maneuvers 
achievable  using  the  ROE  formulas  and  input-shaping  control  profiles,  namely,  re-phasing  or  rendezvous  maneuvers  with 
along-track  control  only. 

The  new  analytical  solutions  in  terms  of  ROEs,  and  particularly  their  simplifications,  as  done  in  the  input-shaping  case,  can  be 
implemented  onboard  small  spacecraft  with  limited  computing  capabilities,  allowing  them  to  autonomously  compute  guidance 
trajectories  of  several  kinds.  In  addition,  optimization  routines  running  on  the  ground  could  be  envisioned,  using  as  initial  guess 
the  here  proposed  formulas.  The  ROEs  analytical  solutions,  simplified  when  applying  the  problem  specific  constraints  (such  as 
fixed  value  for  the  thrust,  duration  of  thrust,  etc.)  can  be  used  as  the  core  for  fast  direct  methods  of  numerical  optimization,  the 
advantage  being  the  a-priori  satisfaction  of  the  constraints  provided  by  the  nature  of  the  solutions. 
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